# Direct outputs eligible for public release

# GLOBAL SETTINGS --------------------------------------------------------------

options(
    scipen = 999,
    digits = 16,
    max.print = .Machine$integer.max,
    show.error.locations = TRUE,
    warn = 1
)

RNGkind("L'Ecuyer-CMRG")
seed <- 818675309L
set.seed(seed) # setting main seed

# PACKAGES ---------------------------------------------------------------------
library(data.table)
library(ggplot2)
library(scales)

# PACKAGE SETTINGS -------------------------------------------------------------

# data.table
setDTthreads(threads = 1L)
options(datatable.print.class = TRUE, datatable.print.keys = TRUE)
# so that printing the data.table also shows the variable type on top

# BEGIN FILE -------------------------------------------------------------------

# Read in all estimates

claim_oasi <-
    readRDS(
        "~/estimation-output/wealth_effect_estimates_has_ssinc_62_to_64.rds" # nolint
    )
claim_oasi <-
    claim_oasi[
        model == "reduced_form" & ref_event_time %in% c(-7:5)
    ]

# Add omitted
temp <- data.table(model = "reduced_form", income_quartile = 5)
temp[, ref_event_time := -2]
temp[, estimate := 0]
temp[, cluster_se := 0]

claim_oasi <-
    rbindlist(list(claim_oasi, temp), use.names = TRUE)
temp <- NULL
rm(temp)

one_yr_exit <-
    readRDS(
        "~/estimation-output/wealth_effect_estimates_nonemployed_62_to_64.rds" # nolint
    )
one_yr_exit <-
    one_yr_exit[
        model == "reduced_form" & ref_event_time %in% c(-7:5)
    ]
one_yr_exit[, outcome_label := "One-Year Exit"]

two_yr_exit <-
    readRDS(
        "~/estimation-output/wealth_effect_estimates_two_years_nonemployed_62_to_64.rds" # nolint
    )
two_yr_exit <-
    two_yr_exit[
        model == "reduced_form" & ref_event_time %in% c(-7:5)
    ]
two_yr_exit[, outcome_label := "Two-Year Exit"]

five_yr_exit <-
    readRDS(
        "~/estimation-output/wealth_effect_estimates_five_years_nonemployed_62_to_64.rds" # nolint
    )
five_yr_exit <-
    five_yr_exit[
        model == "reduced_form" & ref_event_time %in% c(-7:5)
    ]
five_yr_exit[, outcome_label := "Five-Year Exit"]

exits <-
    rbindlist(list(one_yr_exit, two_yr_exit, five_yr_exit), use.names = TRUE)

# Add omitted
temp <-
    CJ(
        outcome_label = unique(exits$outcome_label),
        model = "reduced_form",
        income_quartile = 5
    )
temp[, ref_event_time := -2]
temp[, estimate := 0]
temp[, cluster_se := 0]

exits <- rbindlist(list(exits, temp), use.names = TRUE)
temp <- NULL
rm(temp)

exits[
    ,
    outcome_label :=
        factor(
            outcome_label,
            levels = c("One-Year Exit", "Two-Year Exit", "Five-Year Exit")
        )
]

# Produce Figure B.11a
# Claiming OASI (percentage points)
claim_oasi[, estimate := estimate * 100]
claim_oasi[, cluster_se := cluster_se * 100]

figure_b_11_a <-
    ggplot(
        aes(
            x = ref_event_time,
            y = estimate
        ),
        data = claim_oasi
    ) +
    geom_line() +
    geom_ribbon(
        aes(
            ymin = estimate - (1.64 * cluster_se),
            ymax = estimate + (1.64 * cluster_se)
        ),
        linetype = 0,
        alpha = 0.2,
        show.legend = FALSE
    ) +
    theme_bw(base_size = 13) +
    theme(panel.grid.minor = element_blank()) +
    scale_x_continuous(
        breaks = seq.int(from = -7L, to = 5L, by = 1L),
        expand = expansion(mult = c(0.0025, 0.0025))
    ) +
    scale_y_continuous(breaks = seq.int(from = -3, to = 15, by = 3)) +
    labs(x = "Event Time (ℓ)", y = "Event Study Estimate (pp)") +
    theme(
        legend.title = element_blank(),
        legend.position = c(0.0825, 0.0825),
        legend.background = element_rect(fill = "white", color = "grey"),
        strip.text.x = element_blank(),
        strip.background = element_blank(),
        legend.key = element_rect(NA),
        legend.spacing.y = unit(2, "mm"),
        legend.key.height = unit(4, "mm"),
        legend.key.width = unit(4, "mm"),
        legend.margin = margin(t = 0, r = 0.1, b = 0.15, l = 0.1, unit = "cm")
    ) +
    scale_color_viridis_d() +
    scale_fill_viridis_d() +
    guides(color = guide_legend(ncol = 1, byrow = TRUE)) +
    coord_cartesian(ylim = c(-3, 15))

ggsave(
    plot = figure_b_11_a,
    filename = "~/paper/figures/figure-B.11a.png",
    width = 6,
    height = 4,
    dpi = 600
)

ggsave(
    plot = figure_b_11_a,
    filename = "~/paper/figures/figure-B.11a.tif",
    width = 6,
    height = 4,
    device = "tiff",
    dpi = 600
)

# Housekeeping
claim_oasi <- NULL
figure_b_11_a <- NULL
rm(claim_oasi, figure_b_11_a)

# Produce Figure B.11b
# Labor market exit (percentage points)
exits[, estimate := estimate * 100]
exits[, cluster_se := cluster_se * 100]

figure_b_11_b <-
    ggplot(
        aes(
            x = ref_event_time,
            y = estimate,
            color = factor(outcome_label),
            fill = factor(outcome_label)
        ),
        data = exits
    ) +
    geom_line() +
    geom_ribbon(
        aes(
            ymin = estimate - (1.64 * cluster_se),
            ymax = estimate + (1.64 * cluster_se)
        ),
        linetype = 0,
        alpha = 0.2,
        show.legend = FALSE
    ) +
    theme_bw(base_size = 13) +
    theme(panel.grid.minor = element_blank()) +
    scale_x_continuous(
        breaks = seq.int(from = -7L, to = 5L, by = 1L),
        expand = expansion(mult = c(0.0025, 0.0025))
    ) +
    scale_y_continuous(breaks = seq.int(from = -3, to = 15, by = 3)) +
    labs(x = "Event Time (ℓ)", y = "Event Study Estimate (pp)") +
    theme(
        legend.title = element_blank(),
        legend.position = c(0.117, 0.8825),
        legend.background = element_rect(fill = "white", color = "grey"),
        strip.text.x = element_blank(),
        strip.background = element_blank(),
        legend.key = element_rect(NA),
        legend.spacing.y = unit(2, "mm"),
        legend.key.height = unit(4, "mm"),
        legend.key.width = unit(4, "mm"),
        legend.margin = margin(t = 0, r = 0.1, b = 0.15, l = 0.1, unit = "cm")
    ) +
    scale_color_viridis_d() +
    scale_fill_viridis_d() +
    guides(color = guide_legend(ncol = 1, byrow = TRUE)) +
    coord_cartesian(ylim = c(-3, 15))

ggsave(
    plot = figure_b_11_b,
    filename = "~/paper/figures/figure-B.11b.png",
    width = 6,
    height = 4,
    dpi = 600
)

ggsave(
    plot = figure_b_11_b,
    filename = "~/paper/figures/figure-B.11b.tif",
    width = 6,
    height = 4,
    device = "tiff",
    dpi = 600
)

# Housekeeping
one_yr_exit <- NULL
two_yr_exit <- NULL
five_yr_exit <- NULL
exits <- NULL
figure_b_11_b <- NULL
rm(one_yr_exit, two_yr_exit, five_yr_exit, exits, figure_b_11_b)